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Abstract. We study time evolution of a subsystem's density matrix under 
unitary evolution, generated by a sufficiently complex, say quantum chaotic, 
Hamiltonian, modeled by a random matrix. We exactly calculate all coherences, 
purity and fluctuations. We show numerically that the reduced density matrix 
r^ ■ can be described in terms of a noncentral correlated Wishart ensemble for which 

O ! we are able to perform analytical calculations of the eigenvalue density. Our 

description accounts for a transition from an arbitrary initial state towards a 
random state at large times, enabling us to determine the convergence time after 
which random states are reached. We identify and describe a number of other 
interesting features, like a series of collisions between the largest eigenvalue and 
the bulk, accompanied by a phase transition in its distribution function. 
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1. Introduction 

Isolated systems are usually idealizations that is nevertheless often valid to a good 
approximation due to small coupling to external degrees of freedom. Frequently 
though, one would like to describe the dynamics of a subsystem, even if the coupling 
to "external" degrees of freedom is not small. This is of obvious importance for 
foundations of statistical mechanics where in certain simple situations, for instance in 
equilibrium, one can derive the stationary distribution of a subsystem under rather 
general assumptions. Much more difficult, but also more interesting, is the problem of 
describing the subsystem in a nonequilibrium setting, for instance, its evolution with 
time. In the present work we achieve just that. 

We completely describe how the statistical properties of a reduced density matrix 
evolve in time, provided we have a unitary evolution on the whole system. Because we 
want to avoid system-specific features of a particular Hamiltonian we use a random 
matrix model. The Hamiltonian of the whole system is given as a random matrix, 
in our case from a unitary ensemble, although straightforward generalizations are 
possible also to other symmetry classes. It is well established [T that random matrix 
Hamiltonian's are useful to describe certain statistical properties of quantum chaotic 
systems, meaning that our results should apply to a large class of dynamical systems. 
For such random matrix model of a reduced density matrix we derive exact analytic 
results for all coherences as well as for purity. Calculating spectral properties of 
the reduced density matrix exactly though is still too difficult. With the goal of 
nevertheless obtaining analytical result in the limit of large systems we statistically 
describe the reduced density matrix by a Noncentral Correlated Wishart Ensemble (n- 
CWE) [21 E], i.e., a Wishart Ensemble [4] (WE) with nonzero average matrix elements 
and correlations between rows or columns. Doing this we are able to calculate the 
eigenvalue density of the reduced density matrix, thereby fully describing all invariant 
properties of the reduced density matrix at any time. As an example, calculating the 
average value of the largest eigenvalue we observe a series of its "collisions" with the 
bulk of the spectrum, each time accompanied by a phase transition in its distribution 
function. 

Various Wishart ensembles have been studied in different contexts. The 
uncorrelated Wishart ensemble with a centered distribution, and its fixed-trace 
generalization, has received a lot of attention from physicists. Even a partial list 
of its applications includes diverse fields, like mesoscopic systems [5^, analysis of time 
series in quantitative finance [6j [7] or of EEC signals 8^, high-energy physics [9] 
or communication engineering [101 . It describes also properties of random quantum 
states, important in quantum information theory. A number of interesting phenomena 
has been discovered, like phase transitions in the distribution of extreme eigenvalues or 
of Reny entropies [ill 1121 US] . Various generalizations of Wishart ensemble have also 
been studied T^. On the other hand, n-CWE, being the ensemble that applies to our 
work, though it appeared early in the literature jl5) . is rather unexplored compared 
to WE. It has been mostly of interest in mathematical statistics [21 [T71 [TSl (THj or in 
signal processing [TQ]. This work is the first application of n-CWE to physics, thereby 
opening many new research directions. It turns out that for large matrices correlations 
are small and one deals with an uncorrelated Noncentral Wishart Ensemble (n-WE) or, 
due to universality, in some cases with a Correlated Wishart Ensemble (CWE) which 
has also been of interest in recent years [201 1211 1221 • Because the natural ensemble 
in our setting is the n-CWE this should also give a fresh impetus for more detailed 
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mathematical studies of n-CWE. 

Our results also have implications for the much studied quantum protocols to 
generate random quantum states. An important question is how long does it take to 
generate a random quantum state. For different random circuit protocols |23| . that 
can be described as a Markovian process, it has been proved [24] that the convergence 
time (number of two-qubit operations; one operation per time step) scales with the 
number of qubits n as '--^ n^ (as opposed to previous lower bound [25]). However, 
from an experimental point of view it might often be easier to implement a chaotic 
Hamiltonian than a specific fixed random circuit protocol. A proof of convergence 
time in such setting is still missing, although it has been argued [26] that it should 
have the same scaling as for random circuits. Our general result provides this missing 
link, deriving the convergence time for a large class of protocols based on sufficiently 
complex Hamiltonians. We also show that during evolution the statistical properties 
of states are well described by the n-WE. One interesting observation is that, because 
the convergence to the asymptotic invariant distribution is not monotonic, evolving for 
longer time can, somehow counter-intuitively, worsen the randomness of states. How 
the spectrum of the reduced density matrix changes with time has been numerically 
observed in a random protocol [26J and in a dynamical system [27j . 

The outline is the following: in Section [2] we shall present our model for a reduced 
density matrix and obtain some exact results. In SectionOwe approximate the model 
from the previous section by a Wishart ensemble, obtaining exact results for the 
eigenvalue density in the limit of large system size. In Section [4] we verify theoretical 
results by numerical simulations. In Sec. [5] we discuss the implication of our results 
for the convergence rate towards random states and in Sec. [6] we conclude. 

2. The model and averages over the unitary Haar measure 

Let us first define our random matrix model for a reduced density matrix. Dynamics 
on the whole system is given by a unitary propagator U* = exp {—\Ht) acting on a 
N (S) M dimensional Hilbert space (without loss of generality we assume N < M). We 
are interested in the reduced dynamics on an N dimensional factor space obtained 
by tracing over the "environment" of dimension M . The state at later time |V'(i)) is 
given through the unitary evolution \ip{t)) = f/*|V'(0)), and the reduced density matrix 
is p{t) = tr e|V'(0)(V'(OI- The only assumption we need in our exact derivations in 
this section is that the eigenvectors of H are well described by a random unitary 
matrix so that we will be able to perform averages over unitary Haar measure of these 
eigenvectors, denoted by ()u, while the eigenvalues of H can be arbitrary. We shall 
shortly call such a model a random density matrix model (RDMM). 

As a concrete example of RDMM, by which we will numerically check our 
theory, we will use a Hamiltonian from a Gaussian Unitary Ensemble (GUE) in 
which matrix elements are independent complex Gaussian random numbers [T]. We 
fix normalization, i.e., the variance of random numbers, by (|-ffjfcp) — j^, giving 
the asymptotic (large NM) spectral span of 4, while the Heisenberg time is 27V Af. 
Note that for a specific dynamical system normalization will be of course different, 
however, this results in a trivial rescaling of two relevant time scales: the shortest 
time scale is the inverse spectral span (the largest energy scale in the system), while 
the longest, i.e., the Heisenberg time, is given by the inverse of the mean level spacing 
(the smallest energy scale in the system). The initial state on the total system is 
an arbitrary superposition |V'(0)) — J2i j/^i,!'(0)b)l^)i where \j) are basis states of a 
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central system and {ly) are basis states of an "environment" . The initial density matrix 
isp(o) ^A{0)A1{0). 

We shall now calculate some exact results for the RDMM. Writing [/* ^ = 
Vj^mdiagjexp (— ii?mt)}V^*^ in terms of eigenvectors Vj^k and eigenenergies Em, and 
performing average over unitary Haar-measure distributed eigenvectors V of H, we 
get for the average expansion coefficients A{t) at time t, 

(A(i))u = ft A(0), f.^-L-Y^ exp i-iE,t), (1) 

j 
where ft depends only on the eigenenergies Ej of H. Details of the calculation can 
be found in the appendix. We could also perform averaging over the spectrum, ()e, 
obtaining the average {ft)E = gt, and (|/ip)E = h^- The eigenvector-average of the 
density matrix at time t involves averages over monomials consisting of two matrix 
elements Vj.m and two conjugate T^*,„, resulting in (see the appendix) 

All averages over the Haar measure of Vj,m can be evaluated exactly by the method 
developed in Ref. [28]. Even though the expression ([2]) is only a starting point for 
our calculations it has already important consequences. For instance, it gives all 
coherences (off-diagonal elements of p) which are of interest in studies of decoherence. 
From Eq. ([2]) it is also simple to derive a quantum map p^^' — ?► p{t) for our RDMM, 
which is diagonal in a local basis. This generalizes a recent result [IH] for qubits 
(A'' = 2). Furthermore, Eq.© will be used to define an appropriate Wishart ensemble 
in Section [3l 

Average purity I{t) = {tr E[p^{t)])u can also be evaluated exactly in our setting 
in which the eigenvectors of H are given by a random unitary matrix. The calculation 
involves averages over products of 4 elements Vj^m and 4 conjugate ones, see the 
appendix. The Weingarten function on permutations of 4 elements needed for such 
a calculation [28] can be found in Ref. [29]. After straightforward, but tedious 
calculation, we arrive at the final exact result for the average purity for the RDMM 
and a product initial state (any A^, M > 1 and t) , 

m =I, + B {M'\ft\'+J^vt + |/2tP - 4|/tp) , (3) 

with /r = i^^'li being the average purity of complex random states [33], N = NM, 

vt = [f?]*f2t + fHf2t]* and B = (^|^)7X^|Y)7_;^Li) ■ Purity for the GOE case at certain 
times has been studied in Ref. |34j . Interestingly, for t — > cx) purity does not go 
towards I^ but instead [3S] to Ip -I- j^f/fj'+l^rj^i) - For any finite size and any time, 
even an infinite one, one does not reach perfectly random states! For discrete-time 
evolution described by circular ensembles this has been observed in Ref. |36j . 

In a special case when H is from a GUE, and for large M, we can simplify I{t) 
even more. The infinite size limit of ft averaged over GUE eigenenergy spectrum 
is gt = Ji{2t)/t. In terms of gt we can write in the leading order in A/" that 
h^ = g^ + 0{\/N), (IAHe = gt + 0{\/N^), and {vt)j, = 2g^,g^t + 0{\/N^). Using 
these we obtain the first three leading orders of purity for GUE Hamiltonian 

m = gt + ^^-^'^^^ + "^^ + M931^9tl ^ ^(i/_^.). (4) 

^ ' ^* NM NM ^ y ' J W 

We shall now check our exact result for purity ([3]) against result of numerical 

simulation of RDMM for a GUE Hamiltonian. Average purity is shown in Fig. [TJ 
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where one can observe a perfect agreement between the exact purity ^ and results 
of numerical simulation. 




Figure 1. Purity for a RDMM, H G GUE and a product initial state. Full curve 
is the theory, Eq. Q, circles are numerics; N = M = 256. Note "revivals" of 
purity around t ^ 2.5 and t fa 4.2. 

We have derived an exact expression for (/c(i))u and for the purity in our RDMM. 
The exact calculation of the spectral density on the other hand still seems to be too 
difficult. With that in mind we will make one further simplification. To nevertheless 
be able to obtain analytical results and to give a complete statistical description for 
the subsystem dynamics we shall in the next section describe the RDMM by a simpler 
statistical ensemble, namely by a noncentral correlated Wishart ensemble (n-CWE). 
Such description of the RDMM is only approximate, however, as we shall demonstrate, 
for large dimensions N, M, on which we shall focus in the present work, it becomes 
exact. The ensemble properties are self- averaging in this limit, meaning that ensemble 
averaged quantities describe also properties of one particular instance of p{t). We leave 
small size-cases for future studies. 

3. Wishart ensembles 



The idea is to describe the RDMM in an approximate way by a suitable Wishart 
ensemble. Wishart ensemble parameters are chosen so that the low-order moments 
are equal to those of the RDMM, for instance, equal to the average reduced density 
matrix ([2]). Higher order moments of the RDMM and those of the Wishart ensemble 
will in general differ, however, we expect that the contributions of higher moments 
become negligible for large system sizes. We shall verify this assumption by numerical 
simulations. 
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3.1. Generalized Wishart ensembles as a model for RDMM 

The most general Wishart ensemble is n-CWE. Member matrices W of n-CWE are 
defined as 

W = ZZ\ Z = Y+^X, (5) 

through an auxiliary matrix 2", in which F is a fixed matrix determining an offset 
("nonzero average of matrix elements" in n-CWE), and X is an iV x M matrix of 
independent complex Gaussian variables with zero average and variance {\Xjk\'^) — 
j^ , where () denotes an ensemble average over Gaussian numbers in X. ^ is a positive 
definite symmetric matrix which takes account for the correlations among the rows 
of Z. Note that we have an ensemble average (W) = YY'' + jj£,. We remark here 
that initially CWE has been introduced as a model for the covariance matrices where 
one considers Y = 0. Additionally, statistical independence of the matrix elements is 
imposed for the simplicity from technical point of view. On the other hand, for RDMM 
both conditions are a priory not necessary, or even not valid (i.e., average matrix 
elements of Y are nonzero and there are nonzero correlations) . These generalizations 
therefore make n-CWE more suitable candidate for approximating RDMM. 

We now fix parameters of n-CWE by demanding that Z ^ describes A{t) Q 
while W should describe density matrix p{t) ^. Specifically, we have that {Z) equals 
the average expansion coefficients {A{t)) ([1]), while the covariance matrix ^ is on the 
other hand fixed by equating (W) = cr'^C + lfftPp^°^ = {p{t))v, where ct^ = 1/iV. This 
therefore gives us n-CWE parameters 

Y^gtAiO), a'C = {p{t))v~\9tfp^°^. (6) 

Using Eq.([2]) we can see that for large N and M correlations are small, that is ^ is 
proportional to the identity matrix, (_ — {I ~ h"^)! + 0{\/{NMY). In the rest of the 
paper we shall focus on the case of large N and M and we will in turn approximate 
n-CWE by an uncorrelated n-WE for which ^ ~ 1. 

Even though n-WE is a simpler ensemble than n-CWE compact analytical 
expressions are still difficult to obtain. For instance, the ensemble averaged eigenvalue 
density, (p{x)), is not known in a simple form j30l . Fortunately, in some of the 
cases we discuss below, we can approximate n-WE with that of a simpler CWE. An 
approximation of using CWE to describe n-WE has been used before [31]. CWE 
description of the RDMM is obtained by setting 

r = 0, a2e=(p(t))u. (7) 

Note that in this case of CWE description correlations ^ are not small because they 
must also account for a nonzero p^^' term in p{t) ^ which is in n-CWE approximation 
taken care of by a nonzero Y . For most of the parameters and quantities studied in 
this work we conjecture that the level density is the same for CWE and n-WE, in the 
leading order. In the absence of rigorous analytical treatment, we leave this conjecture 
based on some numerical experiments. A non-rigorous argument for the universality, 
at least for not too small times, is that for the spectral properties of p{t) only its 
moments matter and not for instance expansion coefficients of \'ip{t)), which is fixed 
by a nonzero Y in n-CWE. Similar universality has also been discovered in some other 
random matrix ensembles |16) . 
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3.2. Results from the CWE 

We proceed by repeating analytical results for the eigenvalue density of CWE in the 
limit of large N,M, for more details the reader can consult Ref. [3T]. For a given 
spectrum of £, and for large N, M, the ensemble averaged density of eigenvalues for 
CWE, {p{x)), can be determined in terms of its resolvent G{z), 

{G{z)) = 1 Tr (^z - ^ [« - 1 + z{G{z))] dj , (8) 

where k = M/N, as (p(A)) = ^G{X — ie)/7r, for positive infinitesimal e. For a simple 
spectrum of ^ the Eq. (|5]) can be solved analytically. For instance, for an initial product 
state one has (again, for large sizes) ^jt = Sjt + (1 — ^jk)r, where r = h^. In this case, 
for large r the largest eigenvalue centered at (Ai) and with Gaussian fluctuations [15] 
around, is separated from the bulk X2,...,n [H] (Aj are in decreasing order), and equal 
to 

^^^^^^.(Nr + l-r^iNrn + l^^ (9) 

NrK 
This equation is valid only if the eigenvalue is separated from the bulk, i.e., if 
r > {N^/k)~^. It can also be generalized to a case when more than one eigenvalue is 
separated [21] [16]. The distribution of eigenvalues of p{t) in the bulk is given by the 
Marcenko-Pastur law [32] with a rescaled variance (1 — r)cr^: 



V(A+-A)(A-A.) 
(mp(A))^. 27rA(l-r)a^ ' ^''^^ 

where A± = (1 — 7')(t^(k^^''^ ± l)'^, are the minimum and maximum cut-off of the 
eigenvalue density for the bulk. 

As long as the largest eigenvalue is separated from the bulk we can also 
analytically calculate its variance, though not using the CWE approximation [57] . 
If we start with an initial state in which one eigenvalue of p^°^ is separated from 
the rest - the simplest example is an initial product state - the eigenvector \xt) 
corresponding to this eigenvalue is almost the same as |a;f=o)- Therefore, the variance 
can be approximately calculated from the variance of the matrix element (a;o|/3(i)|a;o)u- 
For an initial product state and RDMM such an average can be evaluated exactly (for 
A^ — > oo) by the same method as used for the exact calculation of the average purity 
(O. The final resuh is 

a (Ai) « a (poo) = j^^ + 0{j^^j^,). (11) 

Note that for product initial states and short times the largest eigenvalue is much 
larger than all the others and will therefore govern properties of the reduced density 
matrix. As we shall see in the next section (e.g.. Fig. [3]), the largest eigenvalue is 
separated from the bulk also at times between collisions of the largest eigenvalue with 
the bulk, that is not necessarily at short times, provided A^ is large enough (one 
therefore has to let A^ — > oo before letting time to infinity). 

4. Numerical comparison 

In this section we shall verify our theoretical results for the eigenvalue density of p{t) 
for the RDMM, obtained from the approximation by the CWE, Eqs. (I10l9p . against 
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Figure 2. (Color online) Density of eigenvalues for RDMM (points) and theory 
(lines), Ea. l|9ll0| l. for different values oft and a product initial state. In the insets 
we show the largest eigenvalue, whose density is Gaussian with theoretical average 
(|9]l and variance mill . 



results of full numerical simulation of RDMM with GUE Hamiltonian. To obtain 
RDMM results we expanded short-time propagator exp {—iHAt) into a power series, 
keeping a sufficient number of terms, and acting with it on pure state. Repeating 
this procedure pure state at a later time is obtained, from which we then calculate 
the reduced density matrix. For largest iV's used the above power-expansion method 
is faster than directly diagonalizing H. We shall also compare these results with the 
numerical simulation of n-WE in order to assess the error made in approximating n- 
WE results with those of OWE. n-WE simulations are performed by straightforwardly 
generating the reduced density matrix using Wishart matrices ^, calculating Y and 
^ according to Eq.([6]), using theoretical values for GUE ensemble of gt = Ji(2i)/i, 
where Ji is the Bessel function, and ht = gt, both holding for large NM (for large 
NM used n-CWE and n-WE results are indistinguishable on the scale of figures). 

Three different cases of initial states will be considered: (i) product initial 
state, (ii) entangled initial state with only two nonzero Schmidt coefficients and (iii) 
entangled initial state with all N nonzero Schmidt coefficients. 



4-.1. Product initial state 

We start with the case (i), product initial state. In Fig. [2] we show the eigenvalue 
density of the RDMM obtained by numerical simulation. These are compared with 
the theory for the bulk density (fTO|) in the main plot and Gaussian distribution with 
the mean given by ([9|) and fluctuations by (fTTj) in the insets. Good agreement can 
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be seen between the theory and RDMM results. We can see that the leading order 
of the CWE theory, giving a rescaled Marcenko-Pastur density in the bulk, agrees 
well with the RDMM simulations, apart from small discrepancy around the edges of 
the spectrum. Theory for n-WE (data not shown) would be on the scale of the plot 
indistinguishable from the CWE theory, except near the bulk edges, where it would 
correctly account for small deviations of RDMM from the Marcenko-Pastur density. 

Next, we proceed by a more detailed comparison of (Ai). For an uncorrelated WE 
and large matrices (Xj) have been calculated in Ref. [M]- In Fig. |3]we show results 
of numerical simulation of RDMM, numerical simulation of n-WE, and analytical 
expression ([9]) obtained from CWE. We can see a nice agreement between all three. 




Figure 3. Average values of two largest eigenvalues for the RDMM simulation 
(points) and n-WE simulation (full curves) for an initial product state. The inset 
shows the "collision" between Ai and the bulk (filled triangles are theory ^ 
obtained for CWE). All is for N = 256, k = 1. 

In Fig. 13] similar comparison is made for the variance of the largest eigenvalue. 

In the plot of (Ai) and (A2), Fig. [31 one can nicely see an interesting "collision" 
between the largest isolated eigenvalue and the bulk, shown by the 2nd largest 
eigenvalue at its edge. Collision does not result in an exact degeneracy (see for instance 
the inset in Fig. [6|) but only in a very close encounter, with (Ai) — (A2) ^ l/N"^ . 
Collisions happen approximately at times when Ji{2t) has zeroes. Few things can be 
observed. First, there are interesting "ripple" effects at the time of collision (see Figs.|3] 
andEl). More interestingly, at the time of the collision there is a phase transition in 
the distribution of the largest eigenvalue [ITl [18] from a Gaussian before the collision, 
to a Tracy- Widom distribution [39 at the collision. This is demonstrated in Fig. [5] 
Opposed to phase transitions in the distribution itself [13] we here have a cascade of 
phase transitions in time. Note that the frequency of these collisions is proportional 
to the spectral span (4 for our CUE normalization), which is the fastest time scale in 
the system. 
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Figure 4. The variance of Ai for RDMM simulation (circles), theory Eq. ljll|l 
(full curve) and a fixed-trace n-WE simulation (dotted curve). All is for N = 
256, K = 1. 
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Figure 5. The distribution of the largest eigenvalue for n-WE (points) and 
theoretical Gaussian/Tracy- Widom distributions (full curves), indicating a phase 
transition with time. All is for A'^ = 256, k = 1. 



4-2. Entangled initial state 

Next, we proceed to show results for two entangled initial states. The first one we 
choose is an entangled state with only two nonzero Schmidt coefficients, |'0(O)) = 
•\/l/4|0)|0) -I- \/3/4|l)|l). The eigenvalue density again agrees perfectly with n-WE 
simulations (data not shown). We only show a plot of three largest eigenvalues, Fig. [6] 
One can again see that n-WE describes our RDMM well and that the largest eigenvalue 
is given by the analytical expression (jH]). 

Finally, we choose a very entangled initial state, in which all Schmidt coefficients 
are nonzero, |^(0)) — ^ y^2j /N{N + l)\j)\j)- The eigenvalue density is shown in 
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Figure 6. Three largest eigenvalues for N = M = 128 and initial entangled 
state with two nonzero eigenvalues Ai = 3/4 and A2 = 1/4. Points are RDMM 
simulation, full curves n-WE theory obtained by numerical simulation and filled 
triangles in the inset Eq.l(9]| with r = ^/ij . 



Figs. [71 [5] respectively for short tinie and for long time, again displaying an agreement 
with the n-WE theory. The oscillatory behaviour in Fig. |5]is due to the small system 
size. 
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Figure 7. (Color online) Density of eigenlevels for N = 32, k = 1 and initial 
state with \j = 2j/N(N + 1) for j = 1, ..., A'^; points are RDMM simulation, lines 
arc n-WE theory. Short time behavior showing transition from a "picket-fence" 
spectrum is shown. In Fig. |8] long time behavior is shown. 
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Figure 8. (Color online) Long time behavior of the same data as shown in Fig.[7l 
At t Ri 1.9, the first zero of gt oc Ji(2t)/t, correlations in the Wishart ensemble 
^ are zero and the density is described by the Marcenko-Pastur form. 



5. Convergence to random states 

Because purity can be used as an indicator of convergence to random states our exact 
results for purity in RDMM enable us to comment on the convergence rate towards 
random states. The Eg.® tells us that I{t) — /^ sa gf. The convergence to a random 
state is reached with 1/N accuracy (which is equal to the scaling of /r with N), 
I{t) — Ir ~ l/iV, when we have gt ~ l/N-^^^. Because gt is the Fourier transformation 
of the spectral density its decay time is given by the inverse of the width of the spectral 
density. If the width is independent of the system size, as is the case for GUE, the 
decay time scales as ~ 1. li gt would be an exponential (Lorentzian eigenlevel density), 
such small value of gt would be reached after time r ~ logiV = n. In a dynamical 
system one can decompose the propagator [/* — exp {—iHt) using e.g. Suzuki- Trotter 
formula into basic two-qubit operations. In a nearest-neighbor coupled chaotic system 
the number of two-qubit operations in a decomposition of a short-time propagator 
U scales with the number of qubits n linearly as ^ n. Decomposition up to time 



demands of order ■ 



two-qubit operations. Our result therefore indicates the 



same scaling with the number of two-qubit gates as in random circuits [24j . 

However, if gt is not an exponential there are two interesting issues: (i) if gt 
has zeroes (as is for instance the case for GUE) one can reach random states already 
after a much smaller time r ^ 1. Note that if the spectral span of eigenenergies is 
bounded one necessarily has zeroes in gt. This is very promising as one could decrease 
convergence time by a factor of n compared to known protocols |24j . (ii) On the 
other hand, for a bounded spectrum one will have a singular behavior at the spectral 
edge, for instance a square-root singularity for GUE, resulting in a very slow algebraic 
decay of the envelope of the oscillating gt. This will result in an exponentially large 
time ^ 2" after which these oscillations die out. A bounded spectrum therefore has 
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advantages and disadvantages. 

It is an interesting open question how these issues are reflected in a concrete 
dynamical system. What is clear from our discussion is that what matters is the 
Fourier transformation of the spectral density, as well as the spectral span, which 
trivially rescales all times given in our work. "Engineering" an appropriate spectral 
density one could speed-up the convergence. 

6. Conclusion 

Using numerical simulations and asymptotic arguments we showed that one can 
describe the statistics of the reduced density matrix during unitary evolution by a 
complex Hamiltonian with a noncentral Wishart ensemble. We have analytically 
calculated purity, the average largest eigenvalue and its variance. Interesting collisions 
between separated eigenvalues and a random-matrix-like bulk occur as time progresses, 
leading to phase transitions in the distribution of these eigenvalues. We discuss 
convergence time for a large class of protocols for generation of random states, also 
suggesting a new way to speed-up these protocols. The results presented open many 
new questions in physics as well as in mathematics. For instance, it would be 
interesting to study generalization to orthogonal or symplectic symmetry. 
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Appendix 

In this Appendix we sketch the steps needed in evaluating averages over unitary Haar 
measure, needed in calculation of average reduced density matrix ^, average purity 
^ and fluctuations of the largest eigenvalue pT|) . They can also be found in previous 
works, for instance, in Refs. [28, 29 . 

Let us write the propagator U*^,.^ = V,^,„Adiag{exp (-i£',„Ai)}V;r^,„A ^ terms 
of eigenvectors Vj^^^ku and eigenenergies Emx- Summations are implied over repeated 
indices. We use notation with double indices, for instance j'/i, where Latin indices 
are for the central subsystem of interest, being of dimension TV, while Greek indices 
enumerate the environmental basis states, being M in number, over which we trace 
in order to obtain the reduced density matrix. Expansion coefficients Aj^^(t) of the 
state at time t are therefore 

AjAt) = C/j.,fe^^fe,40). (A.l) 

To evaluate the average expansion coefficient we therefore need 

{A,At)h = (^..,mAVi;,„A)u e-'^'-*A,,^(0). (A.2) 

Using the well known result for the Haar average of a product of two matrix elements, 

(V,i/,mAVj,^ „ja)u = 'Tf"jv ' (^-3) 
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we immediately get the Eq.([T]). We abbreviated J\f = NM and use a short notation 
for a product of Kronecker delta symbols, Sj^j^'" = Sj^_^^5j^^^^ ■ ■ ■. 

To calculate the average value of the reduced density matrix, obtained as 
p{t) — A{t)A'^{t), we need the average of a product of four matrix elements of V 
(two from each A). They are equal to 



{Va,b 



^a2b2^a3b3^a4b4 



1 



AA2 



1 



^3630464 
aibia2b2 



r04 640363 
01610262 



1 



^03640463 I ^04630364 
01610262 01610262 



,(A.4) 



where we used a short notation Va^ti = VjiVi,kiij.i with a double-index a^ = jiVi and 
hi = fci/ii. Using this formula in the expression for p{t), summing over all repeated 
indices, and finally simplifying the resulting products of Kronecker deltas, one arrives 
at the average reduced density matrix as given in Eq. ([2]). 

For purity, which is quadratic in the reduced density matrix, we have to average 
over product of 8 Vs. For this one needs the Haar average over a product of 8 matrix 
elements. Exact formula is rather lengthy, however, using the result in Ref. [2S|, 
averages over unitary group can be written in the following compact form 

, ...a' 6' ...6' 



(V^iti • ■■Va^bpV*. 



■■ -V* r., 






a\ ...ap bi ...bp 



^^ Wg(AA, 



'\ 



(A.5) 



where the summation is over all permutations a and r of p elements. Weingarten 
function Wg is a rational function in size M and depends only of the cycle shape 
of the permutation crr"^. The whole evaluation of unitary Haar averages therefore 
boils down to knowing values of the Weingarten function. For p = 2 and p — ?> 
they are given in Ref. [28], for p = 4, needed in purity calculation, they have been 
calculated in Ref. [29] . For permutations of p = 4 elements there are 5 different values 
of Weingarten function (i.e., 5 different cycle shapes), while the sum in Eq. (jA.5l) results 
in a sum of (4!)^ — 576 products of Kronecker deltas. A sheer number of terms and 
the fact that one has to sum all these terms over repeated indices, taking into account 
also exponential factors involving energies, makes calculation rather long and tedious. 
After summation and simplification many terms actually result in the same expression 
so that the final results can be compactly written as given in Eq.Q. 
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